Below is an example of 20 hypothetical 1ha restoration plots (10,000m2, 100m*100m) distributed at random at three reefs in the Townsville region (Davies Reef, Big Broadhurst Reef, Little Broadhurst Reef). The code takes the Allen Coral atlas geomorphology data (see here for more details), filters all reef polygons larger than 1ha in size, randomly selects 20 polygons across the three reefs, and constructs a 100*100m plot around the centroid of each polygon:
# data packages
library(tidyverse)
library(units)
library(foreach)
# spatial packages
library(sf)
# mapping packages
library(tmap)
library(leaflet)
# create function to get make polygons for plots
set_restoration_plot_centres <- function(input = NULL, width = NULL, length = NULL) {
# Calculate the coordinates of the rectangular polygon
x <- sf::st_coordinates(input)[1, 1]
y <- sf::st_coordinates(input)[1, 2]
# set parameters
x_min <- x - (width / 2)
x_max <- x + (width / 2)
y_min <- y - (length / 2)
y_max <- y + (length / 2)
polygon <- sf::st_polygon(list(rbind(c(x_min, y_min), c(x_min, y_max), c(x_max, y_max), c(x_max, y_min), c(x_min, y_min)))) |>
sf::st_sfc(crs = 20353)
return(polygon)
}
habitats <- c("Plateau", "Back Reef Slope", "Reef Slope", "Sheltered Reef Slope", "Inner Reef Flat", "Outer Reef Flat", "Reef Crest")
cairns_geomorphic <- st_read("/Users/rof011/coraldynamics/data/Moore-Arlington-20230906220745/Geomorphic-Map/geomorphic.geojson") %>%
sf::st_transform(crs = sf::st_crs("EPSG:20353")) %>%
dplyr::group_by(class) %>%
dplyr::mutate(habitat_id = paste0(
gsub(" ", "_", class), "_",
sprintf(paste0("%0", ceiling(log10(max(1:length(class)))), "d"), 1:length(class)))) %>%
sf::st_make_valid() %>%
# mutate(habitat_id=as.factor(habitat_id)) %>%
dplyr::group_by(habitat_id, class) %>%
summarise(geometry = st_union(geometry)) %>%
mutate(area = round(st_area(geometry))) %>%
filter(class %in% habitats) %>%
drop_na(class)
## Reading layer `Moore Arlington' from data source `/Users/rof011/coraldynamics/data/Moore-Arlington-20230906220745/Geomorphic-Map/geomorphic.geojson' using driver `GeoJSON'
## Simple feature collection with 3768 features and 1 field
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 145.9182 ymin: -16.89574 xmax: 146.2582 ymax: -16.61326
## Geodetic CRS: WGS 84
cairns_plots_centroids <- cairns_geomorphic %>%
filter(area > set_units(10000,"m^2")) %>%
filter(class %in% c("Reef Slope", "Sheltered Reef Slope", "Back Reef Slope")) %>%
drop_na(class) %>%
st_centroid() %>%
st_cast("POINT")
cairns_plot_outputs <- foreach(i=1:nrow(cairns_plots_centroids), .combine="c") %do% {
tmp <- set_restoration_plot_centres(cairns_plots_centroids$geometry[i], width=100, length=100)
tmp
}
Select layers using the layercontrol on the left, use [ ] for full-screen viewing.
set.seed(1)
cairns_plot_outputs <- cairns_plot_outputs |> st_as_sf() |> mutate(id=paste0("Plot_ID_",seq(1,n(),1)))
cairns_plot_outputs2 <- cairns_plot_outputs[sample(nrow(cairns_plot_outputs), 20), ] # subset to 20
# visualise with thematic maps
map <- tm_view() +
tm_tiles("Esri.WorldImagery", group = "<b> [Seascape]</b> satellite map", alpha = 0.5) +
tm_shape(cairns_geomorphic |> mutate(class=as.factor(class)), id="area", name = "<b> [Seascape]</b> habitats") +
tm_borders(col = "black", lwd = 0.2) +
tm_fill("class", name = "Benthic habitats", id="area", palette = c("Plateau" = "lightskyblue1", "Back Reef Slope" = "darkcyan",
"Reef Slope" = "darkseagreen4", "Sheltered Reef Slope" = "darkslategrey",
"Inner Reef Flat" = "darkgoldenrod4", "Outer Reef Flat" = "darkgoldenrod2",
"Reef Crest" = "coral3"), alpha = 0.6) +
tm_shape(cairns_plot_outputs2, id="plot_id", name = "<b> [Restoration]</b> plots") +
tm_fill(col="darkred", alpha=0.6) +
tm_borders(col="black")
map |> tmap::tmap_leaflet() |>
leaflet::addProviderTiles('Esri.WorldImagery', group = "<b> [Seascape]</b> satellite map", options=leaflet::providerTileOptions(maxNativeZoom=18,maxZoom=100)) |>
leaflet::addProviderTiles('Esri.WorldTopoMap', group = "<b> [Seascape]</b> base map", options=leaflet::providerTileOptions(maxNativeZoom=19,maxZoom=100)) |>
leaflet::addLayersControl(position="topleft", overlayGroups=c(
"<b> [Seascape]</b> satellite map", "<b> [Seascape]</b> habitats", "<b> [Restoration]</b> plots"),
options=leaflet::layersControlOptions(collapsed = FALSE)) |>
leaflet.extras::addFullscreenControl(position = "topleft", pseudoFullscreen = FALSE)